function [invvaralpha] = PseudoInverse(varalpha)
	
nstate = size(varalpha,1);

[u_svd,s_svd,v_svd] = svd(varalpha,'econ');
rank_varalpha = rank(varalpha);


% for i=1:1:nstate
%     if i > rank_varalpha
%         s_svd(i,i) = 0;
%     else
%         s_svd(i,i) = 1/s_svd(i,i);
%     end
% end
for i=1:1:nstate
    if s_svd(i,i) < 1e-4
        s_svd(i,i) = 0;
    else
        s_svd(i,i) = 1/s_svd(i,i);
    end
end

invvaralpha = v_svd*s_svd*u_svd';